home *** CD-ROM | disk | FTP | other *** search
/ Complete Linux / Complete Linux.iso / docs / apps / circuits / spice2g6.z / spice2g6 / spice / Fortran / dcdcmp.f < prev    next >
Encoding:
Text File  |  1989-02-03  |  8.8 KB  |  299 lines
  1.       subroutine dcdcmp
  2.       implicit double precision (a-h,o-z)
  3. c
  4. c     this routine swaps rows and columns in the coefficient matrix accor-
  5. c ding to the numerical pivoting and minimum fillin terms.it then performs
  6. c an in-place lu factorization of the coefficient matrix.
  7. c
  8. c spice version 2g.6  sccsid=tabinf 3/15/83
  9.       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
  10.      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
  11.      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
  12.      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
  13.      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
  14.      5   imynl,imvn,lcvn,nsnod,nsmat,nsval,icnod,icmat,icval,
  15.      6   loutpt,lpol,lzer,irswpf,irswpr,icswpf,icswpr,irpt,jcpt,
  16.      7   irowno,jcolno,nttbr,nttar,lvntmp
  17. c spice version 2g.6  sccsid=miscel 3/15/83
  18.       common /miscel/ atime,aprog(3),adate,atitle(10),defl,defw,defad,
  19.      1  defas,rstats(50),iwidth,lwidth,nopage
  20. c spice version 2g.6  sccsid=cirdat 3/15/83
  21.       common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
  22.      1   nut,nlt,nxtrm,ndist,ntlin,ibr,numvs,numalt,numcyc
  23. c spice version 2g.6  sccsid=status 3/15/83
  24.       common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
  25.      1   xmu,sfactr,mode,modedc,icalc,initf,method,iord,maxord,noncon,
  26.      2   iterno,itemno,nosolv,modac,ipiv,ivmflg,ipostp,iscrch,iofile
  27. c spice version 2g.6  sccsid=flags 3/15/83
  28.       common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts,
  29.      1   lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,itl6,igoof,nogo,keof
  30. c spice version 2g.6  sccsid=knstnt 3/15/83
  31.       common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
  32.      1   gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox,
  33.      2   pivtol,pivrel
  34. c spice version 2g.6  sccsid=debug 3/15/83
  35.       common/debug/ idebug(20)
  36. c spice version 2g.6  sccsid=blank 3/15/83
  37.       common /blank/ value(200000)
  38.       integer nodplc(64)
  39.       complex cvalue(32)
  40.       equivalence (value(1),nodplc(1),cvalue(1))
  41. c
  42.       call second(t1)
  43.       if (ipiv.le.0) go to 12
  44.       if (idebug(11).le.0) go to  3
  45.       call dmpmat(6hdcdcmp)
  46.       idebug(11)=idebug(11)-1
  47.     3 do 10 i=2,nstop
  48.       no=0
  49.       loc=nodplc(jcpt+i)
  50.     5 if (loc.eq.0) go to 7
  51.       no=no+1
  52.       loc=nodplc(jcpt+loc)
  53.       go to 5
  54.     7 nodplc(numoff+i)=no
  55.       no=0
  56.       loc=nodplc(irpt+i)
  57.     8 if (loc.eq.0) go to 9
  58.       no=no+1
  59.       loc=nodplc(irpt+loc)
  60.       go to 8
  61.     9 nodplc(nmoffc+i)=no
  62.    10 continue
  63.    12 n=1
  64. c
  65. c     find next pivot
  66. c
  67.    15 n=n+1
  68.       if (ipiv.gt.0) go to 20
  69.       if (idebug(13).le.0) go to 17
  70.       call dmpmat(6hdcdcm2)
  71.       idebug(13)=idebug(13)-1
  72.    17 if (idebug(14).le.0) go to 18
  73.       if (mode.ne.2) go to 18
  74.       call dmpmat(6hdcdcm3)
  75.       idebug(14)=idebug(14)-1
  76.    18 if (idebug(15).le.0.or.icalc.le.10) go to 19
  77.       call dmpmat(6hdcdcm4)
  78.       idebug(15)=idebug(15)-1
  79.    19 nxti=n
  80.       nxtj=n
  81.       go to 120
  82. c
  83. c     search the coresponding column for max entry
  84. c
  85.    20 vmax=0.0d0
  86.       loci=n
  87.    25 loci=nodplc(irpt+loci)
  88.       if (loci.eq.0) go to 50
  89.       i=nodplc(irowno+loci)
  90.       if (i.lt.n) go to 25
  91.    30 if (dabs(value(lvn+loci)).le.vmax) go to 25
  92.       vmax=dabs(value(lvn+loci))
  93.       go to 25
  94.    50 if (vmax.gt.pivtol) go to 60
  95.       write(iofile,51) n,vmax
  96.    51 format('0*error*:  maximum entry in this column at step ',i4,' (',
  97.      1   1pd12.6,') is less than pivtol')
  98.       igoof=1
  99.       return
  100.    60 epsrel=dmax1(pivrel*vmax,pivtol)
  101.       if (n.ge.nstop) go to 200
  102.       if (ipiv.le.0) go to 120
  103. c
  104. c     pivoting on the diagonal
  105. c
  106.       minop=100000
  107.       nxti=0
  108.       do 70 i=n,nstop
  109.       i1=nodplc(irswpf+i)
  110.       j1=nodplc(icswpf+i)
  111.       ispot=indxx(i1,j1)
  112.       if (ispot.eq.1) go to 70
  113.       if (dabs(value(lvn+ispot)).lt.epsrel) go to 70
  114.       nop=(nodplc(numoff+i)-1)*(nodplc(nmoffc+i)-1)
  115.       if (nop.ge.minop) go to 70
  116.       minop=nop
  117.       nxti=i
  118.       nxtj=i
  119.       if (minop.le.0) go to 95
  120.    70 continue
  121.       if (nxti.le.0) go to 75
  122.       if (nxti-n) 120,120,100
  123. c
  124. c     pivoting on the entire matrix
  125. c
  126.    75 do 90 i=n,nstop
  127.       loc=i
  128.    80 loc=nodplc(jcpt+loc)
  129.       if (loc.eq.0) go to 90
  130.       j=nodplc(jcolno+loc)
  131.       if (j.lt.n) go to 80
  132.       if (dabs(value(lvn+loc)).lt.epsrel) go to 80
  133.       nop=(nodplc(numoff+i)-1)*(nodplc(nmoffc+j)-1)
  134.       if (nop.ge.minop) go to 80
  135.       minop=nop
  136.       nxti=i
  137.       nxtj=j
  138.       if (minop.le.0) go to 95
  139.    90 continue
  140.       if (nxti.gt.0) go to 95
  141.       write (iofile,92)
  142.    92 format('0*abort*:  pivot not in dcdcmp')
  143.       igoof=1
  144.       go to 200
  145.    95 if (nxti.eq.n.and.nxtj.eq.n) go to 120
  146.       if (nxti.eq.n) go to 105
  147. c
  148. c     a pivot has been found
  149. c
  150.   100 load=nodplc(irswpf+nxti)
  151.       lr=nodplc(irswpf+n)
  152.       nodplc(irswpf+nxti)=lr
  153.       nodplc(irswpr+lr)=nxti
  154.       nodplc(irswpf+n)=load
  155.       nodplc(irswpr+load)=n
  156.       noff=nodplc(numoff+nxti)
  157.       nodplc(numoff+nxti)=nodplc(numoff+n)
  158.       nodplc(numoff+n)=noff
  159.       if (nxtj.eq.n) go to 110
  160.   105 load=nodplc(icswpf+nxtj)
  161.       lc=nodplc(icswpf+n)
  162.       nodplc(icswpf+nxtj)=lc
  163.       nodplc(icswpr+lc)=nxtj
  164.       nodplc(icswpf+n)=load
  165.       nodplc(icswpr+load)=n
  166.       noff=nodplc(nmoffc+nxtj)
  167.       nodplc(nmoffc+nxtj)=nodplc(nmoffc+n)
  168.       nodplc(nmoffc+n)=noff
  169.   110 call swapij(nxti,n,nxtj,n)
  170.       nxti=n
  171.       nxtj=n
  172. c
  173. c     calculate contribution from nxti, nxtj and find fill-ins
  174. c
  175.   120 if (n.ge.nstop) go to 200
  176.       n1=nodplc(irswpf+nxti)
  177.       n2=nodplc(icswpf+nxtj)
  178.       locnn=indxx(n1,n2)
  179.       if (ipiv.le.0 .and. dabs(value(lvn+locnn)).lt.pivtol) go to 220
  180. c
  181. c     down col j
  182. c
  183.       locr=nodplc(irpt+locnn)
  184.   125 if (locr.eq.0) go to 180
  185.       i=nodplc(irowno+locr)
  186.       value(lvn+locr)=value(lvn+locr)/value(lvn+locnn)
  187.       locc=nodplc(jcpt+locnn)
  188. c
  189. c     for each column element look up row nxti
  190. c
  191.   130 if (locc.eq.0) go to 170
  192.       j=nodplc(jcolno+locc)
  193. c
  194. c     check for fill-in (i,j)
  195. c
  196.       if (ipiv.le.0) go to 135
  197.       call sizmem(jcpt,isize1)
  198.       call reserv(i,j)
  199.       call sizmem(jcpt,isize2)
  200.       if (isize1.eq.isize2) go to 135
  201.       call extmem(lvn,1)
  202.       nttbr=nttbr+1
  203.       value(lvn+nstop+nttbr)=0.0d0
  204. c
  205. c     locate element (i,j)
  206. c
  207.   135 if (j.lt.i) go to 145
  208.       locij=locc
  209.   140 locij=nodplc(irpt+locij)
  210.       if (nodplc(irowno+locij).eq.i) go to 155
  211.       go to 140
  212.   145 locij=locr
  213.   150 locij=nodplc(jcpt+locij)
  214.       if (nodplc(jcolno+locij).eq.j) go to 155
  215.       go to 150
  216.   155 value(lvn+locij)=value(lvn+locij)-
  217.      1                  value(lvn+locc)*value(lvn+locr)
  218.   160 locc=nodplc(jcpt+locc)
  219.       go to 130
  220.   170 locr=nodplc(irpt+locr)
  221.       if (ipiv.le.0) go to 125
  222.       nodplc(numoff+i)=nodplc(numoff+i)-1
  223.       go to 125
  224. c
  225. c     reduce nmoffc for each element in col nxti
  226. c
  227.   180 if (ipiv.le.0) go to 15
  228.       locc=nodplc(jcpt+locnn)
  229.   185 if (locc.eq.0) go to 15
  230.       j=nodplc(jcolno+locc)
  231.       nodplc(nmoffc+j)=nodplc(nmoffc+j)-1
  232.       locc=nodplc(jcpt+locc)
  233.       go to 185
  234. c
  235. c     done
  236. c
  237.   200 if (ipiv.eq.0) go to 210
  238.       if (idebug(17).le.0) go to 202
  239.       call dmpmat(6hdcdcm5)
  240.       idebug(17)=idebug(17)-1
  241.   202 call matloc
  242.       rstats(49)=rstats(49)+1.0d0
  243.       ipiv=0
  244.       if (lvlcod.eq.2) lvlcod=3
  245.       ifill=nttbr-nttar
  246.       perspa=100.0d0*(1.0d0-dble(nttbr)/dble(nstop*nstop))
  247. c
  248. c  calculation of operation count (operation := `*' or `/'):
  249. c
  250. c     noffr := off-diagonal elements in row, not including diagonal,
  251. c                counting only those elements in the remainder matrix
  252. c     noffc := off-diagonal elements in column, not including diagonal,
  253. c                counting only those elements in the remainder matrix
  254. c
  255. c     then we have
  256. c
  257. c        lu decomposition     requires sigma(2,nstop-1) {noffc + noffc*noffr}
  258. c        forward substitution          sigma(2,nstop-1) {noffc + 1}   +   1
  259. c        backward substitution         sigma(2,nstop-1) {noffr}
  260. c
  261. c     which sums to
  262. c
  263. c               sigma(2,nstop-1) {noffc + noffc*noffr + (noffc+1) + noffr} + 1
  264. c         or
  265. c               sigma(2,nstop-1) {noffc*(noffr+2) + noffr + 1}   +   1
  266. c
  267.       iops=1
  268.       nstop1=nstop-1
  269.       do 205 i=2,nstop1
  270.       noffr=nodplc(numoff+i)-1
  271.       noffc=nodplc(nmoffc+i)-1
  272.       iops=iops+noffr+noffc*(noffr+2)+1
  273.   205 continue
  274.       rstats(20)=nstop
  275.       rstats(21)=nttar
  276.       rstats(22)=nttbr
  277.       rstats(23)=ifill
  278.       rstats(24)=0.0d0
  279.       rstats(25)=nttbr
  280.       rstats(26)=iops
  281.       rstats(27)=perspa
  282.       go to 215
  283.   210 if (idebug(18).le.0) go to 212
  284.       call dmpmat(6hdcdcm6)
  285.       idebug(18)=idebug(18)-1
  286.   212 if (idebug(19).le.0.or.icalc.le.10) go to 215
  287.       call dmpmat(6hdcdcm7)
  288.       idebug(19)=idebug(19)-1
  289.   215 call second(t2)
  290.       rstats(50)=rstats(50)+t2-t1
  291.        return
  292.   220 ipiv=1
  293.       write(iofile,221) n,nxti,nxtj,iterno,time
  294.   221 format(' pivot change on fly: n= ',i5,' nxti= ',i5,' nxtj= ',
  295.      1       i5,' iterno= ',i5,' time= ',1pd12.5)
  296.       rstats(49)=rstats(49)+1.0d0
  297.       go to 20
  298.       end
  299.